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Abstract. A layered system of two-dimensional planes containing fermionic polar molecules can potentially 
realize a number of exotic quantum many-body states. Among the predictions, are density-wave instabil- 
ities driven by the anisotropic part of the dipole-dipole interaction in a single layer. However, in typical 
multilayer setups it is reasonable to expect that the onset and properties of a density-wave are modified by 
adjacent layers. Here we show that this is indeed the case. For multiple layers the critical strength for the 
density-wave instability decreases with the number of layers. The effect depends on density and is more 
pronounced in the low density regime. The lowest solution of the instability corresponds to the density 
waves in the different layers being in-phase, whereas higher solutions have one or several adjacent layers 
that are out of phase. The parameter regime needed to explore this instability is within reach of current 
experiments. 

PACS. 03.75.Ss Degenerate Fermi Gases - 05.30.Fk Fermion systems and electron gas - 67.85.-d Ultracold 
gases, trapped gases 



1 Introduction 

After the great successes of cold atomic gas physics us- 
ing neutral atoms with short-range interaction ^[2] , many 
groups have now set their goals on obtaining ultracold 
samples of polar molecules that have an anisotropic long- 
range interaction [3,4,5,6,7,,8,..9,10j. These can, however, 
lead to strong losses and the design of experimental ge- 
ometries that reduce these effects are now becoming a re- 
ality. In particular, the use of two-dimensional geometries 
can reduce losses and at the same time very interesting 
many-body phases in both single- and multilayer configu- 
rations have been proposed [ni[T2l[T8l[T4lira[THl[T7l[T8l[T9l 
[ ^ I ^ [ ^ [25 1 [M 1 [ ^ E5 1 [?7] . One such proposal concerns the 
potential instability of a single two-dimensional layer with 
polar fermions toward the formation of density-waves as 
the polarization of the molecules with respect to the layer 
plane is varied 28,29J. However, the systems of current ex- 
perimental interest are not single-layer [lOj , and the effect 
of adjacent layers is therefore of concern. 

Using linear response within the the random-phase ap- 
proximation, we consider how interlayer interactions in- 
fluence the density-waves instability and how the critical 
strength is modified by interlayer terms. In order to esti- 
mate the effects of exchange terms, we use many-body 
local field factors. This approach has been successfully 
applied to electron systems. We find that the instability 
is enhanced by the presence of in-phase density-waves in 
neighboring layers. The effect depends on the density of 
fermions in each layer and is most pronounced in the low 



density limit where the critical value is inversely propor- 
tional to the number of layers. The latter effect is largely 
insensitive to the inclusion of exchange terms, Fermi sur- 
face deformation, or changes in the effective mass. The 
density-wave instability will therefore occupy a larger re- 
gion of the zero-temperature phase diagram for a multi- 
layered system as compared to a single layer system. 



2 Linear Response and Effective Interaction 

We consider a multilayer system of fermions with dipole 
moment D and mass m confined in planes parallel to the 
xy-plane and separated by the distance d. In the direction 
normal to the planes, all dipoles reside in the lowest quan- 
tum level which we take to be a Gaussian of width w, i.e. 
(j){z) cx exp(— z^/2?x;^). The dipole moments D are aligned 
by an external field forming an angle 9 with respect to 
the normal of the planes and with a projection onto the 
planes which is parallel to the x-axis. The experimental 
setup is illustrated in Fig. [l] Two dipoles separated by r 
interact with the potential V{r) = D^{1 — 3cos^ 9,-d)/r^ 
where 6'rd is the angle between D and r. We assume that 
the layers all have the same density n of fermions. 

To obtain the instabilities of the multilayered system 
we use linear response theory and the random-phase ap- 
proximation (RPA) as was done for the case of a single 
layer in [28ti29j . Within the RPA framework, the density- 
wave instability occurs at the poles of the density-density 
response function. To treat several layers we extend the 
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Fig. 1. The experimental setup and the density waves corre- 
sponding to the lowest (a), the next lowest (b), and the highest 
(c) eigenmodes for the three layer case. The dipole moments 
form the angle 6 with respect to the normal of the planes. Their 
projection onto the planes is parallel to the wave fronts. 



induced potential which we write as 



(2) 



where the matrix 'V{q)ij = Vij{q) contains the interac- 
tion between layers i and j. We approximate the system 
response to that of a non-interacting Fermi gas respond- 
ing to both the external and the induced disturbance. We 
then have 

Sp{q,Lo) = x°{q,^) [4>eAq,^) + V{q)Sp{q,Lo)] , (3) 

where x'^iqj^)ij — ^ijXiiqi'^) is the matrix of response 
functions of the non-interacting system which is of course 
diagonal. Combining Eqs. ([!]) and ([3| we arrive at the 
following matrix equation the response function 



x{q,u;)=[l-x"iq,uj)Viq)] ' x\q.Lo) 



(4) 



In the case of a single layer this equation reduces to the 
standard RPA expression for the density-density response 
function. Here we are interesed in density-wave instabil- 
ities in the static limit cj = and we have to determine 
the singularities of By inversion, we see that these 

occur when 



det[/-x"(q)V(q)] =0, 



(5) 



and this is the equation that we will solve below. 

We assume here that the density in each layer is the 
same, so that the non-interacting response functions are 
all the same, i.e. x°(q) = X°(9)i ^i^d are given by 



0, ^ [ d'i 



d^k /(fc + g)-/(fc) 
(2^^)' ^k+q - ^k 



(6) 



where = h^k^/2m and / is the Fermi distribution. 
In the two-dimensional case of interest here we have the 
explicit expression [SD] 



x°(g) = 



1- 



2kp 

q 



e{q-2kp) - 1 



(7) 



where kp is the Fermi momentum and 0{x) is the Heav- 
iside step-function. For simplicity, we ignore any Fermi 
surface deformation due to the dipolar interaction |28) . 
We will briefly comment on the influence of such effects 
in Sec. 13.51 



RPA to a multilayer (or multicomponent) system. We can 
write a general density fluctuation in response to an ex- 
ternal potential, (f>ex, in momentum (q) and frequency (a;) 
space as 

(5p(q,w) = x(q,a;)0g^(q,w), (1) 

where Sp is a vector quantity containing the disturbances 
in each layer as entries. Likewise, is now in gen- 

eral a matrix of response function with entries Xij(q, w). 
The interactions between the various layers produce an 



2.1 Exchange Corrections 

The RPA analysis above neglects the role of exchange in- 
teractions. In the single-component Fermi system we con- 
sider here, the exchange effect can be significant. As an 
example, we note that for a momentum-independent po- 
tential, the exchange correction would completely cancel 
the direct term in a Hartree-Fock calculation. The effec- 
tive dipolar interaction that we discuss in the next sec- 
tion depends, however, linearly on momentum. The ef- 
fects of exchange can be included via the Hartree-Fock 
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RPA approximation. This unfortunately involves a non- 
local interaction making the resulting numerics somewhat 
involved. We will not pursue such calculations here, but 
rather follow the simpler local field factor approach that 
has been very successful for the electron liquid ^1] . It at- 
tempts to include the intrinsically non-local effects of the 
exchange term through the introduction of an effective 
local 'exchange' potential, in similar spirit to the highly 
successful density-functional method. This approach has 
been applied to two-dimensional double-layer electron sys- 
tems (see for example Ref. [32]) which is a system closely 
analogous to the one studied here. 

In the multilayer setup considered here, we must be 
careful when including exchange corrections in the cor- 
rect places. Since we assume that there is no tunneling 
between the layers, the layer index is effectively a spin co- 
ordinate, and we therefore have no exchange corrections 
for the interlayer interaction. This means that only the 
diagonal terms in Eq. Q have to be modified. We write 
a diagonal entry in the lorni 



l-Vo{q) [1-G(q)]x°(g), 



(8) 



where G{q) is the momentum-dependent local field factor. 
There are various more or less sophisticated ways to calcu- 
late this factor through self-consistent numerical methods 
[5T] . As we are only interested in estimating the effects of 
exchange correlations on the density wave instability, we 
will follow a more intuitive approach originally introduced 
by Hubbard [33 . 

The Pauli principle introduces the so-called 'exchange- 
hole' in Fermi systems. For large q, i.e. short length scales, 
the exchange-hole cancels the direct interaction and G{q) — ; 
1 for q — >■ oo. For g — > 0, i.e. for long distance, the exchange 
effect should not play a role and in turn G{q) -> 0. Be- 
tween these limits, the detailed functional form of G{q) 
of course depends on the particular form of the bare po- 
tential Vo{q). As we are only interested in the qualitative 
effects of the exchange correlations, it is sufficient to use 
the simple function G{q) = ^ta.ii~^{q/s) which interpo- 
lates between the q — and q — >■ oo limits above. Here s is 
the natural scale in the problem at hand; we take s — 2kp. 



2.2 Effective Dipolar Interaction 

The direct dipole-dipole interaction has an intra- and an 
interlayer part in our multilayered setup. The Fourier trans- 
form of the former can be written |34| 



Vo{q) = 



2ttw 



-P2{cose) - ^{e,a)F{qw) 



(9) 



where 9=191 and a is the azimuthal angle between the 
wave vector q = (q^, qy) and the projection of D onto the 
plane which is parallel to the x-axis. P2{x) is the second 
Legendre polynomial, and we have defined the function 
F{x) = J\x{\ - erf(a;/y2)] exp(a;^/2) with erf(a;) the er- 
ror function. To obtain this formula, the z-direction con- 
fining the dipoles in the layers have been integrated out. 



The interesting angular dependence of the intralayer in- 
teraction is contained in the function ^(0, a) = cos^ d — 
SVC? 6 cos^ a. This function provides the anisotropy in mo- 
mentum space which is absent at 6 = when the dipoles 
are oriented perpendicularly to the layer. For w ^ d, the 
interlayer interaction can be written as |17| 



Vi{q) = -2TTD^^{9,a)qe 



— dq 



(10) 



This approximation holds very well for small w and devi- 
ates less than 10% for w = 0.2d. 

As argued in [25112^] . the most unstable direction is 
found at a = 7r/2. This is a configuration where the 
density-wave is perpendicular to the x-axis in order to 
reduce the side-by-side repulsion of the dipoles while op- 
timizing the attraction from the head-to-tail setup, see 
Fig. [l] In this case we have ^ = cos^ 0. The first term 
in Eq. ^ which is constant in momentum space can be 
discarded since we are working with a single-component 
Fermi system As discussed in [5S], the critical 

value in a single-layer has some dependence on 9. Here 
we are interesting in the effects of multiple layers and we 
thus fix 6 at cos^ 9 = 1/3, but our results can be eas- 
ily mapped to a different angle through the substitution 
£)2 cos^ 0. 

With the choices above, the intralayer interaction be- 



Vo{q) = -:P^F{qw), 



3\/2ttw 



(11) 



whereas the interlayer interaction is simply multiplied by 



a factor of ^. For wq <^ I, Eq. (|ll|) reproduces the po- 



tential used for the single layer case in [28','2Sr. As we are 
mostly concerned with the effects of multiple layers, we 
will also assume wq <Si 1 for simplicity, i.e. we assume 
that Vo{q) is linear in q. This linear momentum depen- 
dence of the intralayer potential was used in [2 9) to ar- 
gue that the density-wave instability must occur at some 
dipole strength always. We note that the most unstable 
mode is expected to he at q = 2kF (neglecting the effects 
of Fermi surface deformation). For consistency, we must 
therefore have that 2kpw ^ 1. In terms of the density 
of a single layer, n, this condition reads wy/ldirn ^ 1. 
Thus, either the density must be small or the transverse 
confinement strong. In terms of typical physical scales in 
experiments, we have 



2kpw = 7.1: 



Ifiui V lO^cm 2 



(12) 



Using experimentally relevant values [101 d = 0.5/im and 
assuming w/d = 0.1, we find that 2k pw < 1 for densities 
n < 8 - 10® cm~^ which is fulfilled by current experiments. 



3 Instability Conditions 

We first recapitulate the findings for a single layer in the 
RPA neglecting exchange. The instability equation is 



(13) 
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Assuming that the instabihty occurs first a,t q = 2kF, we 
have x° = and V^ilkp) = -^-KB'^kpl'i. We thus 

have the relation 



(14) 



where is the critical dipole strength [S^. We define a 
dimensionless measure of the strength and density; g := 
2mD'^kp jZf? . We thus have the critical value go = 1. 
To include exchange, we need to make the substitution 
Vb(2/ci.) ^ Fo(2fcj.)[l - G{2kF)\ ^ Vo{2kp)/2 and we ob- 
tain go = 2 instead. To highlight the effects of the multi- 
layer setup, we now proceed to discuss the bi- and trilayer 
cases without the 1 — G{q) factors and defer the discussion 
of exchange corrections to Sec. 3.4 



3.1 The Bilayer 

For the case of two adjacent layers we get the following 
algebraic equation from Eq. ^ 



[l-x'{q)Vo{qf~-[x°{q)VMy 



0. 



(15) 



We note immediately that if we set Vi (q) = we recover 
the usual RPA condition for density instabilities. It is also 
clear at this point that the bilayer will have a smaller 
Dc than the single layer above since the Vi{q) term is 
negative. If one considers the interlayer interaction from 
the point of view of induced interactions this is no surprise 
as such interaction are usually attractive at lowest order. 

At g = 2fci? we can solve the equation above and find 
the lowest critical value for a density-wave instability in a 
bilayer 



1 



9b 



1 



< 1 = 50- 



(16) 



For kpd = 1, we get roughly a 12 percent reduction, 
whereas for a lower density of kpd — 0.5 the difference 
is 27 percent. The other solution to the bilayer equation is 
36 = (1 - e~^^'"^)~^, so that gb < 1 < ijb for aU kpd. Solv- 
ing for the corresponding zero eigenmodes of x(g, 0)^^ 
from Eq. Q, we find 







'1 


Sp2 




1 



for gt and 







■ 1 " 


Sp2 




-1 



for ~gb. (17) 



Here, Spi — 6pi{2kp) is the density fiuctuation in layer i. 
We see that the density waves in the two layers are in- 
phase for the lower solution gb and out of phase for the gb 
solution. Thus, the instability is enhanced by the density 
waves in neighboring layers being in-phase gaining more 
attractive head-to-tail energy and minimizing the side-by- 
side repulsion. Likewise, when the density-waves are out- 
of-phase the instability is suppressed. 



3.2 The Trilayer 

The case of three layers produces the algebraic equation 



'2{x'{qmq)r-{x°{q)V2{q))'] 
~2{x\q)V,{q))\\q)V2{q), 



(18) 



where we have introduced the notation V2{q) for the inter- 
layer potential of the two outer layers that are a distance 
2c? apart. Note that V2{q) differs from Vi{q) by a factor 
of exp(— qd) and we thus expect it to be a much smaller 
quantity than Vi{q) at 2k p. 

In light of the above, we therefore first consider the 
simpler case of V2{q) — 0, i.e. we include only nearest- 
neighbor interactions. This means that the equation for 
the instability factorizes and the condition becomes that 
either (l - x^{q)Va{q)) = or 

[(1 - x''{q)Vo{q)? - 2{x°{q)VM)\ = 0- (19) 

Clearly the latter condition produces a lower critical value 
and we find for the trilayer with V2{q) = that 



1 



9t 



1 



/2e- 



-2kFd 



< 9b, 



(20) 



where the asterisk indicates that we include nearest-neighbor 
interactions only. As expected the trilayer has a reduced 
critical value. A naive guess for the trilayer might be to 
multiply the interlayer strength by a factor of 2 and then 
consider it as a bilayer problem. However, our result demon- 
strates that the enhancement is only by a factor of V2. As 
it turns out the trilayer equation has a rather simple an- 
alytic solution. The roots are 



9t 



1 af+2 

2 

1 a^+2+Va''+8a^ 

2 l-a2 
1 

l-a2 



(21) 



where a = eyip{—2kpd) and we have listed them in order 
of increasing magnitude. The top solution is always less 
than one and decreases with kpd (we denote it by gt in the 
following) whereas the others are always larger than one 
and increase with kpd. The corresponding eigenmodes are 



Spi 




1 








" 1 ' 


5p2 




Va2+8-a 
2 




Va'^+8+a 
2 







Sps 




1 




1 




-1 



(22) 



where layer 2 is the one in the middle, see Fig. [TJ Again, 
the lowest solution corresponds to the density- waves in the 
different layers being in-phase with amplitude now being 
the largest for the layer in the middle. The second eigen- 
mode has the middle layer out of phase and of larger mag- 
nitude than the outer layers. This is the same situation as 
the gb solution for the bilayer, only now the out of phase 
effect is more costly. This is also reflected in the fact that 
this solution is always larger than gb for any kpd, whereas 
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the opposite holds for the lowest solution, i.e. gt < gb- 
The last solution is an eigenmode with density waves of 
the outer layers out of phase and no amplitude change in 
the middle layer (within the RPA). We sketch in Fig. [l] 
the density waves for the three eigenmodes. The physi- 
cally relevant mode for the instability is of course the one 
with the lowest critical value corresponding to the density 
waves in the planes being in-phase. 

In Fig. [2] we plot the (lowest) bilayer and trilayer crit- 
ical values at which the density-wave instability appears 
for kpd < 2. The single-layer critical value, go — 1, is ap- 
proached asymptotically, however, for the range plotted 
the multilayer cases are all below that value by at least 
a few percent. As expected the bilayer is always above 
the trilayer value. For kpd — >■ 0, the critical value un- 
dergoes the largest reduction which is a factor of two for 
the bilayer, while for the trilayer it is a factor of 3 (we 
return to this fact below). We also compare the trilayer 
with and without the interaction of the two outer layers, 
V2{q)- When excluding the term, we see a larger critical 
value, gj", for all kpd than when taking it into account 
in gt- The additional attraction of V2{q) thus reduces the 
critical value as one would expect. 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 



kpd 

Fig. 2. Critical value for the appearance of a density- wave 
instability at 2kp for 6 = cos~^(^) and a = ■k/2 excluding 
exchange effects. The full (black) line is for a bilayer system, 
whereas the trilayer with all interaction is shown as a dashed 
(red) line and with only nearest-neighbor interaction as a dot- 
ted (blue) line. 



3.3 Multiple Layers 

For more than three layers we expect similar behavior as 
seen above, i.e. a critical strength that decreases with de- 
creasing kpd. In the limit kpd — we can in fact find the 
exact solution for gj^ for any number of layers, N. Here the 
matrix in Eq. ([5| simplifies considerable since it has 1 — .g 
in all diagonal and ~g in all the non-diagonal entries. It is 



easy to verify that a vector consisting of ones in every en- 
try is an eigenvector of this matrix with eigenvalue 1 — Ng. 
We thus conclude that the system is unstable towards the 
formation of in-phase density waves in all planes for the 
critical coupling strength g^ = 1/N. The limit of very 
small kpd should therefore approach this simple value. 
This limit is clearly seen for the lowest critical value in 
the bi- and trilayer cases above with the corresponding 
eigenmodes approaching one in all entries. If we take this 
limit by reducing d while keeping kp constant, we see that 
for large N the critical value approaches zero as the layers 
come closer. Here the (non-interacting) system is equiv- 
alent to that of N spins moving in two dimensions. The 
response function Eq. ([6| is then multiplied by a factor of 
N which reduces the critical value by a factor 1/7V. 



3.4 Exchange effects 

The exchange correction has to be included in the diago- 
nal terms of the response function only as discussed in the 
previous section. This means that the eigenvectors corre- 
sponding to the critical couplings are the same irrespec- 
tive of whether the exchange effect is included or not. In 
the large kpd limit, the off-diagonal terms of the interac- 
tion are negligible, and we thus obtain the critical value 
g = 1/{1 — G{2kp)). Without exchange the single-layer re- 
sult is recovered, i.e. g = 1- Using the value G{2kp) = 1/2 
as estimated in the previous section, the limit is a factor 
of two larger. These arguments make it clear that the ef- 
fects of exchange are more pronounced in limit of large 
kpd where the intralayer correlations dominate. However, 
the interlayer correlations dominate for small kpd which 
means that the exchange effects are insignificant in this 
limit. This means that one of our main results, the 1/N 
scaling of the critical coupling strength for kpd <^ 1, still 
holds when exchange is included. 

In Fig. [3] we show numerical solutions for the lowest 
critical values for iV = 2, 3, 10, 20, and 30, when neglect- 
ing (lower full (blue) lines) and including exchange (upper 
dashed (black) hnes). The expected decrease of ^at with 
N is clearly seen both with and without exchange correc- 
tions. For example, at kpd = 1 the critical value is 0.88 
for N = 2, whereas for TV = 30 it is 0.77 when neglecting 
exchange. The numbers are 1.57 for = 2 and 1.23 for 
A^ = 30 when including exchange. This trend continues 
for higher A^. Note that all the limits discussed above are 
clearly confirmed by the numerics. In particular the re- 
sults with and without including the local field factor to 
account for exchange approach each other as kpd becomes 
small and when A^ grows. 

The corresponding eigenmodes all have the density 
waves in the layers in-phase. There are also other solutions 
with larger critical coupling strengths as for the bi- and 
trilayer cases. The eigenmodes can be analyzed in similar 
fashion and one finds that the solutions can be organized 
according to the number of adjacent layers that are out of 
phase with each other with the lowest solution (plotted in 
Fig. [3]) fully in-phase across all layers. We speculate that 
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these higher modes correspond to collective modes in the 
striped phase. This will be examined in the future. 



02 1 
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Including exchange---.^ 
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Neglecting exctiange 

/ 




N=2, 3, 10, 20, and 30 







qI \ \ \ \ \ 1 

0.5 1 1.5 2 2.5 3 



kpd 

Fig. 3. Same as Fig.jijbut for iV = 2, 3, 10, 20, and 30 layers 
from top to bottom. Notice the limit at kpd — > which is 
5JV ^ l/N. 



3.5 Effective mass and Fermi surface deformation 

Finally, we briefly address the question of influence of ef- 
fective mass and Fermi surface deformation caused by the 
dipolar interaction. For the single-layer case, these correc- 
tions have been calculated in Ref. ^8. ; it was found that 
these terms pushes the critical value up by about 25% for 
cos^ = 1/3 and about 18% for = 0. In our setup this 
factor has to be included as a prefactor of in Eq. (|4|, 
i.e. it effectively amounts a redefinition of our g. The ne- 
glect of effective mass and deformation effects means that 
our results represent lower bounds. Note again that the 
eigenvectors for the unstable modes are unaltered by these 
corrections. 



and kpd < 1.5 where the system is not unstable towards 
the formation of density waves. This is valuable for the 
study of other phases like superfluidity which persists to 
smaU U dg. 





N=2, 10, and 30 




Density Wave 




Including exchange 


/ 




Neglecting exchange 





qI \ \ \ \ 1 

0.5 1 1.5 2 2.5 



kpd 

Fig. 4. Phase diagram aX 6 = cos^^(^) and a — 7r/2 as func- 
tion of ?7 = mD^ /h^d and kpd for different number of layers 
N = 2, 10, and 30. The density-wave instability occurs above 
the critical lines of which the full (blue) ones neglect while the 
dashed (black) ones take exchange effects into account. 

The regime of validity of the RPA approach augmented 
by the local field factor when applied to dipolar systems 
can be related to the corresponding situation for the elec- 
tron liquid. In the latter case the RPA is known to provide 
reasonable results in the high density limit while it per- 
forms poorly at low densities where the Coulomb to kinetic 
energy ratio, Tj, becomes large [31]. However, for dipolar 
systems the interaction dominates in the high density limit 
whereas the low density limit is weakly interacting. We 
thus expect the RPA to be accurate for low densities and 
weak dipolar strengths, i.e. when g 1. This is precisely 
the case for the large N limit which is our main interest 
in this work. 



4 Phase Diagram 

In the multilayer setup, the interaction parameter, U — 
mD^/h^d, is a convenient dimensionless measure for the 
strength of interactions in the system. In Fig. (|4| we show 
the zero temperature phase diagram in the {U,Kpd) plane 
for = 2, 10, and 30. The more layers, the earlier one 
expects to enter the density-wave regime as before. We 
also see that one can probe the phase diagram by chang- 
ing either the dipole moment or the density of fermions. 
Changing d is also an option. This is, however, somewhat 
harder as U oc 1/d and the lines of constant Ukpd — 3g/2 
are very similar to the lines shown in Fig. Q. We note 
that the inclusion of the exchange term causes an interest- 
ing plateau of the critical values for large N at kpd ^ 1. 
This implies that there can be a large region with U < 2 



4.1 Competing Phases and Finite Temperature 

The zero-temperature phase diagram for density-wave in- 
stabilities presented above needs to be considered in the 
light of other possible ground-states of the layered dipo- 
lar system. In the case of a single layer and in the weak 
couping limit, a p-wave superfluid state was proposed [TH]. 
Likewise, a region of negative compressibility leading to 
collapse of the system was found [18.28, , although this 
happens outside the parameter regime considered here. 
For several layers, the system can become superfluid with 
the Coopers pairs formed between dipoles residing in dif- 
ferent layers [24l[25l[26ll27) . In the strong-coupling limit, 
a single layer of dipoles can also form a Wigner crys- 
tal with a symmetry which depends on how the dipoles 
are aligned with respect to the plane p!4 l [T5 l [T6 t l35] . The 
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presence of bound states in single and multilayer config- 
uration of both fermionic and bosonic dipoles has also 
been given a great deal of attention recently [20,36,37,38, 
[39 l l4 m i4l]. For strong coupling, chains of multiple dipoles 
in bound states could be the relevant degree of freedom 
in the system and the effective interaction of such con- 
stituents should determine the ground-state, and could be 
very different for odd fermionic chains as opposed to even 
bosonic ones. How the phase diagram of single- and multi- 
layer system at zero temperature maps out is an extremely 
interesting topic for future research. 

At finite temperature one expects the physics to be 
governed by the Berezinskii-Kosterlitz-Thouless transition 
(BKT) In the bilayer case the BKT physics is con- 
tained in the pairing order parameter in the weak-coupling 
limit or in a condensate of bosonic dimers in the strong- 
coupling limit |26| . For multiple layers similar dimerized 
phases are expected that are governed by the BKT tran- 
sition |24| . The universal relation for the critical temper- 
ature scales with the superfluid density As the latter is 
proportional to the total density for strong coupling, the 
low density regimes can be difficult to access. We spec- 
ulate that the interlayer interactions could help stabilize 
the low-temperature phases of the system and in turn eas- 
ier to access experimentally as compared to a single layer. 
Again this is a topic for future research. 



5 Conclusions 

We have considered the density-wave instability of dipo- 
lar fermionic polar molecules confined to a stack of two- 
dimensional layers. As the number of layers increases we 
find a reduction of the critical strength to enter the density- 
wave regime at all densities. The corresponding density 
waves are in-phase in all the planes. In the low density 
limit the critical strength even approaches zero as the 
number of layers grow. 

We thank M. M. Parish for numerous discussion and for pro- 
viding valuable references. 
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